#############################################################
#################       Linear Polynomial        #############
###############      October 10, 2018     ##################
####  Rerun: December 15, 2022 - checks out with the paper 
# coefs and se validated in the replication file table1.R

rm(list=ls())
library(foreign)
library(plyr)
library(readstata13)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(stargazer)


data=read.csv("~/Dropbox/Personal Research 2017/replications/karn_nov16.csv")
names(data)
summary(data$Latitude)
summary(data$Longitude)
summary(data$border1)
summary(data$border2)
summary(data$Slope)


#######################
##### Distances #####

#Distance to Mysore-Bombay Border
rd10.mb=data[which(data$NEAR_DIST_border1<7500),] 

table(rd10.mb$border1)

#Distance to Hyderabad-Bombay Border
rd10.hb=data[which(data$NEAR_DIST_border2<7500),] 

table(rd10.hb$border2)


#Mysore
#outcome-health centers
health.mys=lm(health_binary~border1+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd10.mb) #OLS estimation
summary(health.mys)
health.mys.cl=cluster.vcov(health.mys, rd10.mb$dist_name)
health.mys.se=sqrt(diag(health.mys.cl)) #cluster standard errors

#outcome - paved roads
pucca.mys=lm(pucca_binary~border1+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd10.mb) #OLS estimation
summary(pucca.mys)
pucca.mys.cl=cluster.vcov(pucca.mys, rd10.mb$dist_name)
pucca.mys.se=sqrt(diag(pucca.mys.cl)) #clustered standard errors



#Hyderabad
#outcome-health centers
health.hyd=lm(health_binary~border2+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd10.hb) #ols estimation 
summary(health.hyd)
health.hyd.cl=cluster.vcov(health.hyd, rd10.hb$dist_name)
health.hyd.se=sqrt(diag(health.hyd.cl)) #cluster standard errors

#outcome - paved roads
pucca.hyd=lm(pucca_binary~border2+TOT_POP+
        TOT_SC+TOT_ST+Slope+TerrainRug+Latitude+Longitude, data=rd10.hb) #ols estimation 
summary(pucca.hyd)
pucca.hyd.cl=cluster.vcov(pucca.hyd, rd10.hb$dist_name)
pucca.hyd.se=sqrt(diag(pucca.hyd.cl)) #clustered standard errors





stargazer(health.mys, pucca.mys, health.hyd, pucca.hyd, se=list(health.mys.se, pucca.mys.se, health.hyd.se, pucca.hyd.se), digits=3, 
          omit=c("TOT_POP", "TOT_SC", "TOT_ST", "Slope", "TerrainRug", "Latitude", "Longitude"), 
          dep.var.labels=c("Health Centers", "Paved Roads", "Health Centers", "Paved Roads"), 
          covariate.labels = c("Indirect Rule (Mysore)", "Indirect Rule (Hyderabad)", "Constant"),
          add.lines = list(c("Controls","\\checkmark","\\checkmark","\\checkmark","\\checkmark")), 
          omit.stat = c("rsq", "f", "adj.rsq", "ser"))








